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Introduction 


Photographic  coverage  is  the  simplest  method  of  documenting  the  position  and  motion  of 
objects  in  a  laboratory  environment.  In  many  cases,  it  is  the  last  resort  for  measuring  dimensions, 
positions,  and  motion  during  an  event  when  other  measurement  methods  fail  to  provide  needed  data. 
Photogrammetry  is  the  use  of  photographic  images  for  the  measurement  of  lengths  and  coordinates 
of  an  object.  For  many  years,  three-dimensional  (3-D)  analytical  photogrammetry  has  been  an 
essential  tool  in  surveying  and  map  making  used  for  accurate  measurements  of  land  features.  The 
technique  required  the  use  of  expensive  metric  cameras  to  record  aerial  photographs  and  bulky 
comparators  to  reconstruct  the  3-D  coordinates  of  points  on  a  land  surface.  The  analytical  methods 
used  for  the  reconstruction  were  lengthy,  complex,  and  iterative  in  nature,  requiring  the  user  to 
provide  initial  guesses  of  the  coordinates  being  measured.  Although  highly  accurate,  these 
traditional  methods  are  not  readily  applicable  to  laboratory  measurements  where  stereo-base 
photography  is  commonly  used  for  recording  3-D  shapes. 

In  the  early  1970s,  researchers  (Abdel- Aziz  and  Karara,  1971)  simplified  the  complicated 
analytical  procedure  required  for  the  reconstruction  of  3-D  coordinates  from  planar  two-dimensional 
(2-D)  images,  and  introduced  the  direct  linear  transformation  (DLT)  between  the  image  and  object 
spaces.  The  DLT  reduced  the  complicated  calibration  procedures  which  were  required  to  obtain 
inner  and  outer  camera  parameters,  and  eliminated  the  need  for  expensive  metric  cameras,  thereby 
opening  the  door  for  relatively  inexpensive  nonmetric  cameras  to  be  used  in  3-D  photogrammetry. 
Further,  the  DLT  allowed  the  use  of  convergent  stereogrammetry  which  is  suited  ideally  for 
laboratory  close-range  applications. 

The  simplicity  of  the  DLT  method  in  close-range  photogrammetry  spawned  several 
technologies  in  the  fields  of  biomechanics,  sports  medicine,  orthopedics,  automotive  safety,  robotics 
and  other  fields  where  accurate  measurements  of  the  3-D  coordinates  and  motion  of  an  object  were 
required.  In  one  application  where  the  mechanical  properties  of  aortic  membrane  tissue  were  being 
investigated  (Melvin,  Mohan,  Wineman,  1975),  the  DLT  method  was  used  to  accurately  measure 
the  coordinates  on  the  surface  of  aortic  membrane  as  it  inflated  (Alem,  Melvin,  Holstein,  1978).  In 
another  application  (Schneider  et  al.,  1985),  the  coordinates  of  landmarks  on  the  bodies  of  volunteer 
car  drivers  were  measured  in  order  to  describe  the  size  and  shape  of  an  "average"  driver  and  to 
develop  a  crash  manikin  which  is  representative  of  the  50th  percentile  male  driving  population. 

In  a  recent  study,  researchers  at  the  U.S.  Army  Aeromedical  Research  Laboratory 
(USAARL)  investigated  the  effects  of  rear-surface  signature  of  50-caliber  body  armor  as  it  defeated 
a  round.  At  the  same  time,  the  viscoelastic  injury  criterion  was  being  advanced  (Viano  and  Lau, 
1988)  and  had  gained  wide  acceptance.  Briefly,  the  criterion  predicts  an  injury  risk  from  the  product 
V*C  of  two  measured  parameters:  the  rate  or  velocity  (V)  of  deformation  of  the  chest,  and  the 
amount  of  compression  (C)  of  the  chest  in  the  direction  of  impact.  Clearly,  in  order  to  apply  the 
viscoelastic  injury  criterion  or  evaluate  the  effects  of  the  rear-surface  on  injury  production,  it  would 
be  necessary  to  measure  the  time-history  of  the  rear  surface  as  it  deforms.  This  measurement  then 
may  be  used  to  extract  the  speed  of  compression  (V)  of  the  chest  wall,  and  the  amount  of 
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compression  (C)  of  the  chest.  Because  of  the  configuration  of  impact  and  the  speed  of  deformation, 
noncontact  3-D  high  speed  photogrammetry  is  the  only  practical  method  for  accurate  determination 

of  the  V*C  product. 

This  report  describes  the  method  required  to  implement  the  DLT  procedure  for  obtaining 
accurate  static  measurement  of  three-dimensional  coordinates  and  provides  a  computer  program 
which  may  be  used  for  the  calibration  cameras  pointed  to  a  3-D  object  space,  and  for  the 
reconstruction  of  3-D  coordinates  of  point  targets  once  the  field  has  been  calibrated.  The  method 
may  be  extended  to  many  static  and  dynamic  applications,  including  the  measurement  of  the  motion 
of  rear-surface  of  the  .50-caliber  body  armor  as  it  defeats  the  penetration  of  a  round. 


Objectives 

The  objectives  of  this  study  are:  (1)  describe  the  direct  linear  transformation  for  close  range 
photogrammetry;  (2)  present  a  procedure  for  the  calibration  of  a  3-D  object  space  using  the  DLT 
equations,  and  provide  a  computer  program  to  implement  the  calibration  procedure;  and  (3)  present 
a  procedure  for  the  reconstruction  of  coordinates  of  a  point  target  in  the  calibrated  3-D  space  using 
the  DLT  equations,  and  provide  a  computer  program  which  implements  the  reconstruction 

procedure. 


Methods 

The  direct  linear  transformation  treats  the  entire  chain  of  imaging  components  as  a  single 
system  that  reduces  the  3-D  object  space  to  a  2-D  image  plane.  Thus,  we  make  no  distinction 
between  the  camera  proper  and  the  projector.  Instead,  the  entire  camera-projector  system  is 
calibrated  as  a  single  unit.  Further,  because  of  the  use  of  a  single  projector  to  produce  images 
captured  with  different  cameras,  we  refer  to  an  imaging  system  as  a  "camera,"  with  the 
understanding  that  this  "camera"  includes  the  projector  as  well. 

For  any  such  imaging  system  or  "camera,"  it  may  be  shown  that  the  object  space  coordmates 
(x,  y,  z)  of  a  point  target  are  transformed  to  the  image  plane  coordinates  ( u ,  v)  according  to  the  two 

linear  equations: 

u  +  xCj  +  yC2  +  zC3  +  C4  +  uxC9  +  uyCi0  +  uzCn  =  0  ^ 

v  +  xC5  +  yC6  *  zCn  +  Cg  +  vxC9  *  vyCl0  +  vzCn  =  0 


where  (C,,  C2, ...,  Cn)  are  constant  coefficients  which  depend  on  the  imaging  system  being  used  to 
convert  the  3-D  object  (x,  y,  z )  to  a  2-D  image  (u,  v)  coordinates.  The  imaging  system  is  said  to  be 
calibrated  when  its  1 1  DLT  coefficients  are  determined. 
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In  order  to  calibrate  a  camera,  i.e.,  compute  its  1 1  DLT  coefficients,  a  sufficient  number  of 
control  points  are  placed  in  the  field-of-view  of  the  camera  to  produce  at  least  1 1  equations  in  the 
1 1  unknowns.  A  control  point  is  a  target  whose  (x,  y,  z )  coordinates  are  known  precisely  in  the  3-D 
object  space.  When  its  image  is  captured  and  projected  by  a  camera  system,  the  ( u ,  v)  image 
coordinates  of  a  control  point  also  become  known.  Therefore,  the  only  unknowns  in  the  2  linear 
equations  are  11  DLT  coefficients.  Since  2  equations  may  be  written  for  each  control  point,  a 
minimum  of  6  control  points  are  required  to  produce  the  system  of  linear  equations  in  1 1  unknowns. 

The  only  constraint  on  the  placement  of  the  six  control  points  is  that  they  must  not  be  co- 
planar  to  ensures  the  existence  of  a  solution  to  the  system  of  equations.  With  6  control  targets,  12 
equations  in  1 1  unknowns  will  be  written.  Although  one  equation  may  be  dropped  before  solving 
for  the  unknowns,  a  least  square  solution  is  preferred  to  minimize  the  error.  This  implies  that  more 
control  targets  should  be  used  to  produce  an  overdetermined  system  of  equations  in  the  1 1  unknowns 
which  may  then  be  solved  by  least-squares  methods.  Such  practice  reduces  potential  errors  in 
determination  of  ft,  y,  z)  and  the  reading  of  (u,  v)  coordinates,  and  produces  more  accurate  results. 


Given  n  control  point  targets  (n;>6),  we  may  write  2n  equations  in  11  unknowns.  The 
resulting  system  of  equations  is  written  in  matrix  form  as: 


Ti  1 

0  0  0  0 


y,  z<  1 

oooo 


x„  y„  Zn  1 

oooo 


oooo 
*,  yt  zi  1 

oooo 

*,  y,  zt  1 

oooo 

x  y  z  1 

n  *  n  n 


«1*J  “l*  U!Z. 

“l^l  «1Z> 

uix,  u,y,  uiz> 

u,y,  uizi 

UnXn  U»y*  UnZn 

UnXl  Un  yn  UnZn 


1 


(2) 


Equation  (2)  is  an  overdetermined  system  of  linear  algebraic  equations  in  1 1  unknowns 
which  may  be  solved  numerically  with  least  squares  methods.  Figure  1  provides  a  flowchart  of  the 
procedure  which  may  be  followed  for  the  calibration  of  two  cameras. 


Three-dimensional  reconstruction  simply  means  determination  of  the  (x,  y,  z)  position  of  a 
point  in  the  calibrated  3-D  object  space,  given  the  calibration  constants  of  each  camera,  and  the  (u, 
v)  image  coordinates  in  each  camera.  This  is  the  ultimate  product  of  3-D  close  range 
photogrammetry,  that  is,  measures  the  coordinates  of  a  point  target  in  the  laboratory  3-D  space.  To 
accomplish  this,  two  or  more  cameras  are  used  to  capture  the  image  of  the  point  of  interest  at  an 
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instant  of  time.  Once  the  cameras  have  been  calibrated  as  described  earlier,  the  image  of  the  same 
point  target  is  digitized  for  each  camera.  This  generates  as  many  pairs  of  (uh  vk)  coordinates  of  the 
same  target  as  there  are  cameras.  Thus,  for  the  k-th  camera,  equations  (1)  are  written  as: 

Uk  +  X^l,k  +  y^2,k  +  2  C3  +  C4k  +  UkX^9,k  +  U  kyC  10, k  +  UkZ^ll,k  =  ® 

Vk  +  X^5,k  +  yC 6,k  +  Z^7,k  +  ^6,k  +  VkX  ^9,k  +  ^  10,*  +  VkZ^U,k  =  ® 


where  all  terms  are  known  except  the  three  coordinates  (x,  y,  z)  of  the  target  point.  Therefore,  a 
minimum  of  two  cameras  are  required  to  produce  sufficient  number  of  equations  (exactly  four) 
which  may  be  solved  for  the  three  unknowns.  For  practical  reasons,  several  cameras  are  employed 
to  cover  the  same  field-of-view.  Assuming  that  n  cameras  (n^2)  are  available  to  observe  the  same 
point,  equations  (3)  may  be  written  for  each  camera  then  combined  into  the  following  matrix  form: 


(C.,  ♦  ^C91) 

(C*,  +  «,Cuu> 

(Cw  *  Uxcux) 

-  (C4.,  +  Ml) 

(C5jl  *  vjC9 ,) 

(C6,l  +  VlC10.l) 

(Cv  *  v,Cul) 

-  (CM  *  V,) 

(Cu  ♦  «*Cv) 

(C5.*  +  V*C9,*) 

(^2,*  +  Uk^  10,*) 

*  vtC10i) 

(C3,*  +  «*Cn.*) 

(C7>*  ♦  vtCllt) 

X 

<  y  ’  «  - 

Z 

-  (CM  *  Uk) 

-  +  V*> 

(^l,m  +  Um  ) 

(  *  Um  Clo.m ) 

-  (C4>m  *  um) 

.  (C5„  *  vC9„) 

(C6,„  +  V^1M) 

(C7„  +  v.CnJ  . 

.  - (CM *  vj 

Equation  (4)  is  an  overdetermined  system  of  2n  linear  equations  with  (x,  y,  z)  as  the  unknowns,  to 
be  solved  using  least-squares  methods.  Figure  2  provides  a  flowchart  of  the  procedure  which  may 
be  followed  for  the  reconstruction  of  3-D  coordinates  from  the  images  of  two  cameras. 


Results 


Both  calibration  and  reconstruction  phases  involve  the  solution  of  overdetermined  system 
of  linear  equations  using  least-squares  methods.  The  subroutine  LLSQ  listed  in  Appendix  A  is  an 
implementation  of  such  a  method.  It  is  coded  for  the  Microsoft  Fortran  5.0  compiler,  but  easily  may 
be  converted  to  other  programming  languages.  The  stability  of  the  procedure  depends  on  the 
selection  of  the  control  points  and  its  accuracy  on  the  number  of  points  and  the  precision  with  which 
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their  (x,  y,  z)  and  (u,  v)  coordinates  were  determined.  A  necessary  constraint  for  the  stability  of  the 
solution  is  for  the  control  points  to  be  nonplanar. 

The  calibration  procedure  outlined  in  Figure  1  is  implemented  in  the  subroutine  DLTCAL. 
Appendix  B  provides  a  listing  of  the  DLTCAL  subroutine  which  sets  up  the  equations  and  calls  the 
LLSQ  procedure  to  compute  the  1 1  DLT  coefficients.  The  subroutine  is  coded  in  free-form  for  the 
Microsoft  Fortran  5.0  compiler,  but  easily  may  be  modified  for  other  languages.  Input  to  the 
DLTCAL  routine  are  the  number  of  control  points  (at  least  six),  the  laboratory  (x,  y,  z )  coordinates 
of  the  control  targets,  and  the  (u,  v)  image  coordinates  of  each  target  obtained  by  digitizing  the  image 
from  the  camera  being  calibrated.  The  subroutine  returns  the  1 1  coefficients  for  that  camera. 

The  DLTCAL  routine  must  be  called  for  each  camera  being  calibrated.  Since  all  cameras 
must  calibrate  the  same  3-D  object  space,  only  one  array  of  control  (x,  y,  z)  coordinates  must  be  used 
as  input  to  the  calibration  routine,  whereas  the  (u,  v)  array  will  vary  from  camera  to  camera, 
depending  on  the  visibility  of  the  target  from  that  camera  perspective. 

Appendix  C  provides  a  listing  of  subroutine  DLTREC  which  implements  the  3-D  reconstruc¬ 
tion  procedure  outlined  in  Figure  2.  The  DLTREC  subroutine  sets  up  the  equations  and  calls  the 
LLSQ  procedure  to  compute  the  three  unknown  coordinates.  The  routine  is  coded  in  free-form  for 
the  Microsoft  Fortran  5.0  compiler,  but  easily  may  be  modified  for  other  languages.  The  reconstruc¬ 
tion  process,  i.e.,  calls  to  the  DLTREC  routine,  must  be  repeated  for  each  target  point.  Input  to 
DLTREC  are  the  number  of  cameras  (at  least  2),  the  1 1  DLT  coefficients  for  each  camera,  and  the 
image  (w,  v)  coordinates  of  the  target  seen  from  each  camera.  The  routine  returns  the  reconstructed 
laboratory  (x,  y,  z )  coordinates  of  the  target. 


Discussion 


The  DLT  method  described  here  is  a  simple  and  accurate  method  for  measurement  of  3-D 
positions.  It  may  be  applied  to  a  single  point,  to  several  points  of  a  surface,  or  across  several  instants 
of  time  to  perform  3-D  motion  analysis. 

In  practice,  several  cameras  are  calibrated  simultaneously.  Each  camera  produces  its  own 
set  of  DLT  coefficients  which  are  valid  only  for  the  physical  setup  of  the  imaging  system  that 
produced  the  («,  v)  coordinates  employed  in  that  calibration.  Once  the  imaging  system  which 
defines  a  camera  has  been  disturbed  because  of  repositioning  of  the  camera  mount,  refocusing  of  the 
camera  or  projector,  adjustment  of  projection  angle  or  scale,  or  change  of  reference  system  for  the 
image  plane,  the  1 1  DLT  coefficients  become  invalid  and  must  be  determined  anew  by  following 
the  calibration  process.  To  avoid  repeated  calibrations,  cameras  usually  are  locked  into  position  for 
the  duration  of  a  study,  and  are  removed  only  when  the  study  is  completed.  Alternately,  a  fixed- 
based  stereo  system,  with  two  or  more  cameras  attached  to  a  rigid  but  moveable  structure,  may  be 
calibrated  only  once  and  will  produce  3-D  coordinates  relative  to  the  attachment  structure. 
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To  improve  the  precision  of  calibration  and  subsequent  3-D  reconstruction,  and  to  reduce  the 
effects  of  lens  distortion,  the  ( u ,  v)  coordinates  should  be  measured  in  an  image  reference  frame 
whose  origin  is  at  or  near  the  optical  center  of  the  image.  Since  many  digitizing  tablets  report 
coordinates  relative  to  the  lower  left  comer  of  the  digitizing  surface,  it  is  necessary  to  transform  the 
raw  digitized  data  using  a  simple  origin  translation.  This  implies  the  need  for  a  repeatable  procedure 
to  ensure  precise  redetermination  of  the  same  image  center  which  was  used  in  the  calibration.  In  the 
absence  of  fiducial  marks  available  in  metric  cameras,  the  procedure  must  rely  on  hardware  features, 
such  as  the  edges  or  comers  of  the  imaging  gate,  which  are  always  produced  during  the  imaging 
process. 

Researchers  (Abdel- Aziz  and  Karara,  1974)  have  proposed  various  models  to  correct  imaging 
distortions  and  to  account  for  second-order  terms  neglected  in  the  linearization  of  the  object-image 
transformation.  However,  experience  has  shown  that  additional  accuracy  gained  by  the  introduction 
of  second-order  corrections  is  insignificant  when  compared  to  potential  improvement  resulting  from 
precision  and  care  taken  in  the  measurements  of  image  coordinates,  and  does  not  justify  the  increase 
in  equations  and  solution  complexities. 


Conclusions 

A  simple  and  accurate  method  has  been  provided  to  measure  the  3-D  coordinates  of  a  point 
target  in  the  laboratory  using  photographic  coverage  and  the  direct  linear  transformation.  The  full 
procedure  leading  to  3-D  motion  analysis  will  be  presented  in  a  separate  report. 
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Figure  1 .  Flowchart  of  camera  calibration  procedure. 


Camera  1  DLT  Coeff.  I  , - (  3-D  Test  Object  I - N  I  Camera  2  DLT  Coeff. 
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Figure  2.  Flowchart  of  3-D  coordinates  of  a  single  point. 


Appendix  A. 

Listing  of  a  subroutine  to  solve  system  of  linear  equations. 


nnnnwnHnBHW»MHMHww*»»w»nww»*wn«n»»i»i«»i«»i*«ii*i»«iiH"i«iiii!i!i»!t»iiiii*iiiiiiiii!i.ii!iiiiii«iiiiii»i»iiii»iiii"" 

"  Coded  for  MS  Fortran  5.0,  free- form,  by  N.  Alem,  USAARL,  October  1994 

n 


"  Purpose: 

n 

n 

n 

it 

n 

"  Parameters: 

n 

n 


n 

H 

n 

n 

n 

n 

n 

n 

n 

it 

n 

n 

n 

n 

if 

n 

it 

ti 


it 

it 


it 

itnniitfififiiifiiHiiiiHit 


To  solve  the  system  of  linear  equations: 

[A]  *  {X}  =  {B} 

by  minimizing  the  eucledian  norm  of  {b}-[A]*{x},  where  [A] 
is  an  H  by  N  matrix  with  M  not  less  than  N. 


A  ...  M  by  N  coefficient  matrix 

B  ...  M  by  L  right  hand  side  matrix 

M  . . .  number  of  rows  in  matrices  [A]  and  {b} 

N  ...  number  of  unknown,  i.e.,  number  of  columns  of  [A], 

which  is  also  the  number  of  rows  in  matrix  {x} 

L  ...  number  of  columns  in  {x}  and  {b},  (usually  =  1) 

X  ...  N  by  L  solution  matrix  (usually  N  by  1  vector) 

IPIV  ....  Integer  vector  supplied  by  the  use  to  receive 

information  on  column  interchanges  in  matrix  [A] 

EPS  . . .  Input  parameter  which  specifies  relative  tolerance 
for  determination  of  rank  of  matrix  [A] 

IER  ...  Error  parameter.  IER=0  means  successful  solution. 

AUX  ...  Auxiliary  storage  array  of  dimension  max(2*N,L) 

On  return,  first  D  locations  of  AUX  contain  the 
resulting  least  squares. 


SUBROUTINE 

DIMENSION 


LLSQ  (A,  B,  M,  N,  L,  X,  IPIV,  EPS,  IER,  AUX) 
A { * ) ,  B  (*) ,  X{*),  IPIV(*),  AUX (*) 


IF  (M  .LT.  N)  THEN  1  UNDER -DETERMINED  SYSTEM 

IERR  =  -2  l  ERROR  CODE 

RETURN 
END  IF 
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PIV  =0.  i  INITIAL  S(K)  STORE  IN  AUX(K) 

I END  =  0 

DO  K  =  1,  N 

IPIV  (K)  =  K 
H  =  0. 

1ST  =  I END  +  1 
I END  =  I END  +  M 

DO  I  =  1ST,  I END 

H  =  H  +  A  (I)  *  A  (I) 

END  DO 

AUX  (K)  =  H 
IF  (H  .GT.  PIV)  THEN 
PIV  =  H 
KPIV  =  K 
END  IF 

END  DO 

IF  (PIV  .LE.  0.)  THEN 
IER  =  -1 
RETURN 
ENDIF 

SIG  =  SQRT  (PIV)  1  TOLERANCE  FOR  CHECKING  RANK  OF  [A] 

TOL  =  SIG  *  ABS  (EPS) 


1  ZERO -MATRIX  [A] 
I  ERROR  CODE 


n  n  ii  n  n  n  n  ii  ii  n  n  n  n  ii  n  ii  n  n  n  H  n  n  n  n  n  n  n  DECOMPOSITION  LOOP 


n  n  n  n  n  n  n  li  n  n  ll  n  ii  n  n  n  n  ii  n  ll  il  ii  li  ll  li  ll  li  ll 


LM  =  L  *  M 
1ST  =  -M 


DO  K  =  1,  N 

1ST  =  1ST  +  M  +  1 
I END  =  1ST  +  M  -  K 
I  =  KPIV  -  K 

IF  (KPIV  .GT.  K)  THEN  t  EXCHANGE  COLUMNS  K  AND  KPIV 

H  =  AUX  (K) 

AUX  (K)  =  AUX  (KPIV) 

AUX  (KPIV)  =  H 
ID  =  I  *  M 
DO  I  =  1ST,  I END 
J  =  I  +  ID 
H  =  A  (I) 

A  (I)  =  A  (J) 

A  (J)  =  H 
END  DO 
END  IF 
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IF  (K  .GT.  1)  THEN 
SIG  =  0. 

DO  I  =  1ST,  I END 

SIG  =  SIG  +  A  (I)  *  A  (I) 

END  DO 

SIG  =  SQRT  (SIG) 

IF  (SIG  .LE.  TOL)  THEN 
IER  =  K  -  1 
RETURN 
END  IF 
END  IF 

H  =  A  (1ST) 

IF  (H  .LT.  0)  SIG  =  -SIG 

IPIV  (KPIV)  =  IPIV  (K) 

IPIV  (K)  =  KPIV 
BETA  =  H  +  SIG 
A  (1ST)  =  BETA 
BETA  =  1.  /  (SIG  *  BETA) 

J  =  N  +  K 
AUX  (J)  =  -SIG 

IF  (K  .LT.  N)  THEN  1  TRANSFORMATION  OF  MATRIX  [A] 

PIV  =  0. 

ID  *  0 
JST  =  K  +  1 
KPIV  =  JST 

DO  J  =  JST,  N 
ID  =  ID  +  M 
H  =  0. 

DO  I  =  1ST,  I END 
II  =  I  +  ID 
H  =  H  +  A  (I)  *  A  (II) 

END  DO 

H  =  BETA  *  H 
DO  I  =  1ST,  I END 
II  =  I  +  ID 

A  (II)  =  A  (II)  -  A  (I)  *  H 
END  DO 

II  =  1ST  +  ID 
H  =  AUX ( J)  -  A (II )  *  A (II) 

AUX ( J)  =  H 
IF  (H  .GT.  PIV)  THEN 
PIV  =  H 
KPIV  =  J 
END  IF 
END  DO 

ENDIF 


1  SAVE  INTERCHANGE  INFORMATION 
1  COMPUTE  BETA  AND  VECTOR  UK 


1  COMPUTATION  OF  PARAMETER  SIG 


I  TEST  FOR  SINGULARITY 
I  ERROR  (K-l) :  RANK  [A]  LESS  THAN  N 
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DO  J  =  K,  LM,  M  1  TRANSFORM  RIGHT  HAND  SIDE  MATRIX  B 

H  =  0. 

I END  =  J  +  M  -  K 

II  =  1ST 

DO  I  =  J,  IEND 

H  =  H  +  A  (II)  *  B  (I) 

II  =  II  +  1 
ENDDO 

H  =  BETA  *  H 

II  =  1ST 

DO  I  =  J,  IEND 

B  (I)  =  B  (I)  -  A  (II)  *  H 
II  =  II  +  1 
END  DO 
END  DO 

END  DO  l  END  OF  DECOMPOSITION  LOOP 


I  =  N  1  BACK  SUBSTITUTION  AND  INTERCHANGE 

LN  =  L  *  N 

PIV  =  1.  /  AUX  (2  *  N) 

DO  K  =  N,  LN,  N 

X  (K)  =  PIV  *  B  (I) 

I  =  I  +  M 
END  DO 

IF  (N  .GT.  1)  THEN 

JST  =  (N  -  1)  *  M  +  N 
DO  J  =  2,  N 

JST  =  JST  -  M  -  1 
K  =  N  +  N  +  1  -  J 
PIV  =  1.  /  AUX  <K) 

KST  =  K  -  N 
ID  =  IPIV  (KST)  -  KST 
1ST  =  2  -  J 
DO  K  =  1,  L 
H  =  B  (KST) 

1ST  =  1ST  +  N 
IEND  =  1ST  +  J  -  2 
II  =  JST 

DO  I  =  1ST,  IEND 
II  =  II  +  M 
H  =  H  -  A  (II)  *  X  (I) 

END  DO 

I  =  1ST  -  1 

II  =  I  +  ID 

X  (I)  =  X  (II) 

X  (II)  =  PIV  *  H 
KST  =  KST  +  M 
END  DO 
END  DO 
ENDIF 


15 


1ST  =  N  +  1 
TEND  =  0 


1  COMPUTATION  OF  LEAST  SQUARES 


DO  J  =  1,  L 

I END  =  I END  +  M 
H  =  0. 

IF  (M  .GT.  N)  THEN 
DO  I  =  1ST,  I END 

H  =  H  +  B  (I)  *  B  (I) 
ENDDO 

1ST  =  1ST  +  M 
ENDIF 

AUX(J)  =  H 
ENDDO 

IER  =  0 

RETURN 

END 
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Appendix  B. 


Listing  of  subroutii  :  to  calibrate  a  camera. 


infniinnnHniinnnnniiniiniiifi«nifNifi«nnnnnnnnnniiiiniiHiinHniiniinnniinnNNniiiinNnniiifiiiiiiiinnnniiniiii 


"  Coded  for  MS  Fortran  5.0,  free-form,  by  N.  Alem,  USAARL,  October  1994 

n 


"  Purpose: 

n 

n 

n 

n 

"  Input : 

If 

If 

II 

II 

If 

II 

n 

it 

it 

"  Output : 

If 

"  Requires: 

If 

n 

"  Note  1 . 

n 

n 

n 

n 

n 

"  Note  2 . 

H 

n 

n 

n 

"  Note  3 . 

n 

M 

n 


Compute  the  11  Direct  Linear  Transformation  constants  COF 
for  a  camera  system.  The  COF  are  coefficient  of  the  linear 
transformation  which  converts  the  three  corrdinates  (X,Y,Z) 
of  a  point  to  its  image  coordinates  (U,V)  in  that  camera. 


NPT 


. . .  number  of  control  points  (minimum  6) 


XYZ ( 3 ,  * ) 


array  containing  the  3  laboratory  coordinates 
of  the  NPT  points,  known  with  precision. 


UV (2 ,  *) 


COF (11) 


array  containing  the  2  image  coordinates  of  the 
NPT  control  points,  digitized  and  referred  to 
the  optical  center  of  the  image. 

eleven  DLT.  coefficients  for  the  camera. 


Subroutine  LLSQ  to  solve  the  overdetermined  system  of  linear 
equations  with  least  squares  method. 


An  error  flag  is  returned  via  the  *  parameter  to  indicate 
that  the  linear  least  squares  procedure  (LLSQ)  failed  to 
converge  to  a  solution.  In  this  case,  the  parameter  EPS 
which  is  required  inside  the  LLSQ  but  defined  here  should  be 
increased. 


An  error  indication  alBO  occurs  if  the  number  of  points  NPT 
is  less  than  6  or  greater  than  50.  If  more  than  50  points 
are  used  in  the  calibration,  the  dimension  of  AA()  array 
should  be  increased  to  (22  *  NPT),  and  BB ()  to  (2  *  NPT) 

This  calibration  routine  must  be  called  for  each  camera, 
and  must  be  repeated  if  any  of  the  cameras  used  in  the 
system  has  been  altered  or  moved. 


nnnniiHniiiiiiiinnniinniiiiiiiiii 


n 


n 


SUBROUTINE 

REAL *4 

DIMENSION 
REAL  *4 


DLTCAL  (NPT,  XYZ,  UV,  COF,  *) 

XYZ (3 , *) ,  UV (2 , *) ,  COF (*) 

AA(1100) ,  BB (100) ,  AUX (22),  IPIV(ll) 
EPS  /l.E-15/ 
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NEQ  =  2  *  NPT 
12  =  0 

DO  N  =  1,  NPT 


11 

=  N 

12 

=  11  ■ 

+  1 

AA 

(11 

+ 

0 

* 

NEQ) 

— 

XYZ  (1,  N) 

AA 

(XI 

+ 

1 

* 

NEQ) 

= 

XYZ  (2,  N) 

AA 

(11 

+ 

2 

* 

NEQ) 

= 

XYZ  (3,  N) 

AA 

(11 

+ 

3 

* 

NEQ) 

= 

1.0 

AA 

(11 

+ 

4 

* 

NEQ) 

= 

0.0 

AA 

(11 

+ 

5 

* 

NEQ) 

= 

0.0 

AA 

(11 

+ 

6 

* 

NEQ) 

= 

0.0 

AA 

(11 

+ 

7 

* 

NEQ) 

= 

0.0 

AA 

(11 

+ 

8 

* 

NEQ) 

= 

UV  (1,  N)  * 

XYX 

(1, 

N) 

AA 

(11 

+ 

9 

* 

NEQ) 

= 

UV  (1,  N)  * 

XYZ 

(2, 

N) 

AA 

(11 

+10 

* 

NEQ) 

= 

UV  (1,  N)  * 

XYZ 

(3, 

N) 

BB 

(ID 

=  - 

UV  (1 

,  N) 

AA 

(12 

+ 

0 

* 

NEQ) 

s 

0.0 

AA 

(12 

+ 

1 

* 

NEQ) 

= 

0.0 

AA 

(12 

+ 

2 

•k 

NEQ) 

= 

0.0 

AA 

(12 

+ 

3 

* 

NEQ) 

- 

0.0 

AA 

(12 

+ 

4 

* 

NEQ) 

= 

XYZ  (1,  N) 

AA 

(12 

+ 

5 

* 

NEQ) 

= 

XYZ  (2,  N) 

AA 

(12 

+ 

6 

* 

NEQ) 

= 

XYZ  (3,  N) 

AA 

(12 

+ 

7 

* 

NEQ) 

= 

1.0 

AA 

(12 

+ 

8 

* 

NEQ) 

= 

UV  (2,  N)  * 

XYZ 

(1, 

N) 

AA 

(12 

+ 

9 

* 

NEQ) 

= 

UV  (2,  N)  * 

XYZ 

(2, 

N) 

AA 

(12 

+10 

* 

NEQ) 

= 

UV  (2,  N)  * 

XYZ 

(3, 

N) 

BB 

(12) 

=  - 

UV  (2 

,  N) 

END  DO 

CALL  LLSQ  (AA,  BB,  NEQ,  11,  1,  COF,  IPIV,  EPS,  IER,  AUX) 

IF  (IER  .NE.  0)  RETURN  1 

RETURN 

END 
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Appendix  C . 

Listing  of  subroutine  to  reconstruct  a  point  in  3D  space. 


it  ii  n  n  n  n  n  n  n  n  n  n  11  n  11  n  n  n  n  n  n  n  n  n  n  n  n  h  n  n  n  n  n  . .  n  11  n  h  n  n  n  n  "  «  "  •«>*  "  "  »  "  '•  »  "  11  "  »  11  "  11  "  11  "  "  11  "  "  "  "  "  "  *  " 

by  N.  Alem,  USAARL,  October  1994 


Coded  for  MS  Fortran  5.0,  free -form. 
Purpose : 


Compute  the  3  laboratory  coordinates  (x,y,z)  of  a  target  point 
from  its  image  coordinates  (u,v)  pairs  obtained  from  two  or 
more  cameras,  given  11  DLT  coefficients  of  each  of  the  cameras. 


Input : 


NCAM  ...  number  of  cameras  (minimum  2) 

COF (11, *)  ...  eleven  DLT  coefficients  for  each  of  the  cameras. 

UV(2,*)  ...  array  containing  the  image  coordinates  (u,v) 

pairs  of  the  same  point  being  reconstructed, 
observed  and  digitized  in  each  camera  image, 
and  referred  to  the  optical  center  of  the  image. 


Output : 


XYZ (3 ) 


3  laboratory  coordinates  (x,y,z)  of  the  point 
being  reconstructed. 


"  Requires  s 

II 

II 

II 

II 

n 

■  Note  1 . 

it 

n 

n 

ii 

n 

■  Note  2 . 

it 

ii 

IT 


Subroutine  LLSQ  to  solve  the  overdetermined  system  of  linear 
equations  with  least  squares  method.  LLSQ  routine  is  listed 
in  Appendix  A  of  this  report. 


An  error  flag  is  returned  via  the  *  parameter  to  indicate 
that  the  linear  least  squares  procedure  (LLSQ)  failed  to 
converge  to  a  solution.  In  this  case,  the  parameter  EPS 
which  is  required  inside  the  LLSQ  but  defined  here  should  be 
increased. 

An  error  indication  also  occurs  if  the  number  of  cameras  NCAM 
is  less  than  2  or  greater  than  9.  If  more  than  9  cameras  were 
used  to  digitize  the  same  point,  then  the  dimension  of  AA() 
should  be  increased  to  (6  *  NCAM),  and  BB ( )  to  (2  *  NCAM). 


ii  ■■  ii  ii  n  ii  ii  ■■  ii  ii  n  n  ii  ii  n  n  ii  ii  n  ii  ii  H  "  I'  "  "  "  "  »  "  "  "  "  «  n  "  »  ■  "  "  "  n  "  ii  n  ii  n  n  n  11  n  n  n  ■■  »»»•»»'»«»""  "  "  "  »  "  "  "  "  "  "  " 


SUBROUTINE 

REAL* 4 

DIMENSION 
REAL *4 


RECXYZ  (NCAM,  COF,  UV,  XYZ,  *) 

COF ( * ) ,  UV (2 , *) ,  XYZ ( 3 , * ) 

AA (54 ) ,  BB (18) ,  AUX{6),  IPIV(3) 
EPS  /l.E-15/ 
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IF  (NCAM  .LT.  2)  RETURN  1 
IF  (NCAM  .GT.  9)  RETURN  1 

NEQ  =  2  *  NCAM 
12  =  0 

DO  K  =  1,  NCAM 


11  =  K 

12  =  II  +  1 


AA 

(11 

+ 

0 

* 

NEQ) 

= 

COF 

(1, 

K) 

+ 

UV 

(1, 

K) 

* 

COF 

(  9, 

K) 

AA 

(11 

+ 

1 

* 

NEQ) 

= 

COF 

(2, 

K) 

+ 

UV 

(1, 

K) 

* 

COF 

(10, 

K) 

AA 

(11 

+ 

2 

* 

NEQ) 

= 

COF 

(3, 

K) 

+ 

UV 

(1, 

K) 

* 

COF 

(11, 

K) 

AA 

(12 

+ 

0 

* 

NEQ) 

COF 

(5, 

K) 

+ 

UV 

(2, 

K) 

* 

COF 

(  9, 

K) 

AA 

(12 

+ 

1 

* 

NEQ) 

= 

COF 

(6, 

K) 

+ 

UV 

(2, 

K) 

* 

COF 

(10, 

K) 

AA 

(12 

+ 

2 

* 

NEQ) 

= 

COF 

(7, 

K) 

+ 

UV 

(2, 

K) 

* 

COF 

(11, 

K) 

BB 

(il) 

= 

_ 

(  COF 

(4 

,  K) 

- 

UV 

(1, 

K) 

) 

BB 

(12) 

= 

:  - 

(  COF 

(4 

,  K) 

- 

UV 

(2, 

K) 

) 

END  DO 

CALL  LLSQ  (AA,  BB,  NEQ,  11,  1,  XYZ,  IPIV,  EPS,  IER,  AUX) 

if  (IER  .NE.  0)  RETURN  1 

RETURN 

END 
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